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We study quasiparticles of d-wave superconductors in the vortex lattice by self-consistently 
solving the Bogoliubov-de Gennes equations. It is found for a pure d x 2_ y 2 state that: (i) low-energy 
quasiparticle bands in the magnetic Brillouin zone have rather large dispersion even in low magnetic 
fields, indicating absense of bound states for an isolated vortex; (ii) in finite fields with A;f£o small, 
the calculated tunneling conductance at the vortex core shows a double-peak structure near zero bias, 
as qualitatively consistent with the STM experiment by Maggio-Aprile et al. [Phys. Rev. Lett. 75 
(1995) 2754]. We also find that mixing of a d xy - or an s-wave component, if any, develops gradually 
without transitions as the field is increased, having little effect on the tunneling spectra. 



Quasiparticles in superconductors have been a mat- 
ter of a central theoretical interest in condensed matter 
physics, underlying almost all phenomena of supercon- 
ductors. With a growing number of experiments in zero 
magnetic field which support the d x i -^2-wave scenario 
for the high-T c cuprates , efforts have also been made 
to clarify their properties in finite magnetic fields. Ex- 
perimentally, the scanning tunneling microscope (STM) 
is one of the most powerful tools to study them. It was 
Hess et al. who clarified its potential to probe the quasi- 
particles Choosing a conventional s-wave material 
NbSe2, they obtained a beautiful vortex-core spectrum 
with a characteristic zero-bias peak. Its origin was soon 
attributed, through numerical calculations based on the 
Bogoliubov-de Gennes (BdG) equations for an isolated 
vortex |(||7j] , to the bound quasiparticles around the vor- 
tex core predicted long ago by Caroli et al. [|| . Although 
a direct observation of the discrete levels is yet to be per- 
formed, the experiment confirmed the existence of bound 
states in the s-wave pairing and stimulated further the- 
oretical studies on quasiparticles in the vortex states as 
well as similar experiments on high-T c materials. 

The first microscopic theoretical study on a d-wave 
pairing was carried out by Wang and MacDonald us- 
ing a lattice model in finite magnetic fields With 
parameters chosen appropriate for the cuprates, the cal- 
culated differential tunneling conductance at the vortex 
core shows a single broad peak near zero bias, whereas 
that of the corresponding s-wave model exhibits a clearly 
resolved bound-state structure. On the other hand, the 
measurement on YBa2Cu307_a by Maggio-Aprile et al. 
reveals a double-peak structure around zero bias, which 
they attributed by using the s-wave result to a couple 
of bound states split widely due to the large energy gap 
lipf . This discrepancy between the theory and the ex- 
periment raised a couple of questions: (i) whether or not 
bound states exist for the d x 2_ y 2 pairing where the gap 
closes along lines; (ii) why the double-peak structure ap- 
pears. Himeda et al. focused on the strong-correlation 
effect to answer the second question pd| . Choosing the 
two-dimensional t-J model, using a Gutzwiller wavefunc- 
tion, and adopting the Gutzwiller approximation, they 



reported that a double-peak structure may show up in the 
low-doping region due to the induced s-wave component 
which brings a fully opened gap in the core region. How- 
ever, the validity of using the Gutzwiller approximation is 
not well understood here. On the other hand, Franz and 
Tesanovic investigated an isolated vortex of a continuum 
model through a direct diagonalization of a large matrix 
HH]. They have shown that: (i) no bound states exist in 
the d x 2_ y 2 model; (iia) the pure d x 2_ y 2 model also yields 
a single broad peak near zero bias as inconsistent with 
the experiment; (iib) the observed double-peak structure 
may be explained by a d x 2_ y 2 + id xy model where the 
gap opens for an arbitrary direction so that there exist 
bound states in the core region. However, the model (iib) 
suffers from an inconsistency with the experiments that 
the gap is open everywhere even in zero magnetic field. 

With these backgrounds, this letter focuses on a more 
detailed study on a continuum d-wave model. Due to 
technical difficulties, investigations on this fundamen- 
tal model started rather recently, and only the special 
case of an isolated vortex has been considered by now 
p2| p"4| . We here consider the whole range H c \ < H < H C 2 
in a unified way towards a thorough understanding of 
the basic d-wave model, which will also be useful when 
treating lattice models with strong correlations or dis- 
orders in the vortex lattice. The technique we use may 
be called "Landau-level expansion method," developed 
firstly for the s-wave vortex lattice of k ^> 1 |l5|-|l7| , and 
then suitably generalized to include anisotropic pairings 
and the magnetic-field variation |l8| . We thereby clarify 
those effects caused by the formation of the vortex lat- 
tice and seek alternative mechanisms for the double-peak 
structure. We also consider the possibility of mixing of 
a d xy - and an s-wave component into the pure d x 2_ y 2 
state to see how it develops with the field strength and 
how it affects the physical properties. This topic has at- 
tracted much attention recently due to an observation of 
a plateau in thermal conductivity k(H) |jl9), followed by 
a theoretical proposal that it may imply a phase transi- 
tion between the d x 2_ y 2 and the d x 2_ y 2+id xy states [ pfj[ . 
The origin of the plateau is still controversial , and it 
will be worth clarifying theoretically whether or not the 
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transition is really possible. 

Our starting point towards those purposes is the BdG 
equations for the eigenfunctions u s and v* labeled by the 
quantum number s with a positive eigenvalue E s : 



dr' 



U(r,r') A(r,r') 
-A*(r,r') -2f(r,r' 
u s (r) 
-v s *(r) 



u s (r') 
-v*(r') 



(1) 



Here A is the pair potential and 7£ denotes the normal- 
state Hamiltonian giving a quadratic dispersion in zero 
magnetic field |22) ; both are 2x2 matrices to include the 
spin degrees of freedom. The pair potential is determined 
self-consistently by 



A(r,r') = V(r-r')J?(r,r'), 



(2) 



where V denotes the interaction and <& is the order pa- 
rameter defined by 



J?(r,r') 



^ u 8 (r)vj(r')-v 8 (r)uj(r') , E s 



• tanh- 



2k B T 



(3) 



with T denoting the transpose. 

To solve these equations efficiently for the vortex- 
lattice states, we use the Landau-level expansion method 
mentioned above, which is sketched as follows Con- 
sider specifically that the field is along the c axis of a 
layered two dimensional material. We fix the mean flux 
density B rather than the external field H , and expand 
u s and v s in the eigenfunctions of the magnetic trans- 
lation group whose unit cell covers the area <fio/B (0 O : 
the flux quantum). Those eigenfunctions {ipNka{Y)} are 
labeled by the Landau-level index N, the magnetic Bloch 
vector k, and the quantum number a (=1,2) which sig- 
nifies two-fold degeneracy of the orbital states. Now, Eq. 
(|I|) can be solved separately for each (k, a) due to the 
translational symmetry of the vortex lattice; the corre- 
sponding eigenstate is labeled explicitly by s — (v'kaa) 
with v (er) the band (spin) index. On the other hand, Eq. 
(H) is handled via a double expansion of A and $ with re- 
spect to the center-of-mass and the relative coordinates. 
Convenient bases to this end are given by {^Nqi^T - )} 

and {V4!rL( r — r ')}i respectively, with q another mag- 
netic Bloch vector and m the angular momentum along 
the c axis (m = -N, -N+ 1, -N+2, ■ ■ •). The overlap 

between V , ^ q (^)V , Ar r ) „ l ( r - r ') and -0JVka(r)V ; JV'k'a'( r ') 
is calculated rather easily, vanishing unless q = k + k' 
and ct\ = a2- Also worth noting here is that p3| : 

(i) a single q suffices to describe the conventional vor- 
tex lattice due to its broken translational symmetry; 

(ii) No's for the conventional hexagonal (square) lat- 
tice are multiples of 6 (4). Finally (Nm\V\N'm!) is 
estimated most conveniently in the momentum repre- 
sentation. To this end, we expand the pair scatter- 
ing amplitude in zero field: V pp > = (27r) 2 (p|y|p') as 



^pp' = E„ im ' t/(m ' m ' ) (P^')e lm ^ imV -', where <p p de- 
notes the angle of the wavevector p with the a axis. 
Then the matrix element in finite fields is obtained 
as (Nm\V\N'm') = i m ~ m ' (-l) N+N 'V( m ' m ')(p,p')/4Trl% 
with Ib the magnetic length and p=^/N /Is- Once V pp > 
is given explicitly, we can thereby calculate the properties 
of the vortex-lattice state over £f cl <H< H c2 . 
The pairing interaction we use is given by 

N V PP < = -2 [{p/k F ) 2 (p' /k F ) 2 (g x 2_ y 2 cos 2<p p cos 2tp p , 
+g xy sin2<p p sm2<f p <) + g s ] Q(p)Q(p') , (4) 

where Nq is the density of states for both spin in zero 
field, <?'s are the coupling constants, kp is the Fermi 
wavenumber, and O(p) denotes a smooth cut-off [ p"7| : 
Q(p) = e~^ p _fc p)/Pc] with p c a cut-off wavenumber. We 

then utilize the quantity T$ =T$ (gj,p c ) {j = x 2 -y 2 , 
xy, s), defined as the transition temperature in zero field 
for the state j in the absence of other two channels, as 
parameters to measure the strength of the pairing inter- 
action for the channel j. This procedure may be justified 
for the weak-coupling model. Indeed, we find little dif- 
ference in the calculated results among several choices of 
(gj,p c ). Another important parameter in the model is 
the coherence length £o = v f / Ao , where vf is the Fermi 
velocity and Ao denotes the maximum of the energy gap 
at H = T = 0. We will clarify below the fcpCo depen- 
dence of the differential tunneling conductance which is 
proportional to 

N(E,r) = -J2i\u s (r)\ 2 f(E-E s ) 

S 

+ W s (r)\ 2 f(E+E s )}, (5) 

with f{E) = (e £ / feBT + l)" 1 . Although they can be in- 
cluded easily, we here neglect: (a) the spatial variation 
of the magnetic field since k is quite large for the high- 
T c materials; (b) the Pauli paramagnetism for simplicity. 
The square vortex lattice is assumed throughout , and 
the expansion over iV c (the center-of-mass Landau-level 
index) is cut at some A,! nax , typically N™ ax = 8—16, with 
the convergence checked by changing N™^. 

We now present the results on the pure d x 2_ y 2 model 
{g s = g X y = 0). Figure 1(a) shows the quasiparticle spec- 
trum in the magnetic Brillouin zone with /cf£o = 5 for a 
low-temperature and low-field state: T/T c — B j H C 2 — 0.1. 
For comparison, the corresponding one of the s-wave 
model (hexagonal lattice) is presented in Fig. 1(b). We 
can clearly observe the following difference between the 
two. In Fig. 1(b) (s-wave), a qualitative change in the 
dispersion curves is seen around E s — Ao: the curves 
for E s > A are closely spaced with a little dispersion, 
whereas the lower-energy bands are flat and occur in pairs 
with the level spacing of the order of Aq/ef (ep: Fermi 
energy). As already pointed out by Norman et al. [ p7[ , 
these corresponds to the bound core states of an isolated 
vortex with little tunneling probability between adjacent 
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FIG. 1. Quasiparticle dispersion in the magnetic Brillouin 
zone for (a) the d x 2_ y 2 wave (square lattice) and (b) the 
s wave (hexagonal lattice). Parameters are T/T c — 0.1, 
B/H C 2 = 0.1, and fc F £o = 5 

cores; one can check that the number of states for the 
lowest pair of bands is exactly the same as the number of 
vortices. In contrast, the low-energy e?-wave bands in Fig. 
1(a) are densely packed with large dispersion, indicating 
extended nature of the corresponding wavefunctions (this 
feature also shows up in the case of the hexagonal lattice). 
From this comparison, we may conclude that no bound 
states exist for the d-wave model even in the zero-field 
limit of an isolated vortex, in agreement with the result 
of Franz and Tesanovic |l2| . A remark is necessary at the 
same time that, strictly speaking, all the quasiparticles 
in both models become extended in finite fields, as long 
as the vortices are regularly spaced. Thus the difference 
of the low-lying states between the two models in finite 
fields is of quantitative (not qualitative) nature. 

We turn our attention to the differential tunneling 
conductance of the pure d-wave model and study both 
B / H C 2 and &f£o dependence of its spectral line shape. 
Figure 2 shows the quantity at the vortex core calcu- 
lated for T/T c = 0.1 and fc F £ = 5. At the lower field of 
B/H C 2 = 0.05, we see a single broad peak around zero 
bias as previous calculations |j 12 1. It indicates the ex- 
istence of virtually bound quasiparticles in the core re- 
gion. As the field is increased, however, the line shape 
starts changing to show a clear double-peak structure at 
B / H C 2 = 0.3, which is qualitatively consistent with the 
spectrum observed by Maggio-Aprile et al. Q . Though 
not presented here, such a change cannot be seen in the 
spectra obtained for /cf£o = 10 where the single broad 
peak remains essentially unchanged over a wide range of 
the field strength. Hence the splitting may be consid- 
ered as inherent to systems with small values of kp^Q. It 
is explained here in connection with the increase of the 
quasiparticle hopping probability between cores as fol- 
lows. The equilibrium thermodynamic potential at T = 
can be written in terms of Eq. (j|) as 



n 



dr dujLuN{r,u) + Tr<P*A, 



(0) 




FIG. 2. The tunneling conductance of the d a! 2_ !/ 2-wave 



model at the core center for B/H C 2 
T/T c = 0.1 and /fc P £ = 5. 



0.05 and 0.3 with 



where Tr implies integration and summation over the 
space and the spin coordinates, respectively. The sec- 
ond term on the right-hand side is expected not to de- 
pend strongly on the details of the quasiparticle spec- 
trum. The first term, on the other hand, is quite sen- 
sitive to the structure of the local density of states. In 
particular, high density of states at u> — is energeti- 
cally unfavorable so that the system tends to reduce it 
by using any perturbation available. We identify this 
perturbation with the hopping probability of the low- 
energy quasiparticles, whose magnitude may be mea- 
sured by the dimensionlcss distance between adjacent 
vortices, kpls^ ^f^o/ \J B j H C 2- Hence the identification 
tells us that the double-peak structure becomes clearer as 
&f£o (B/H c2 ) decreases (increases), which is consistent 
with our numerical results. Our calculations performed 
for other values of &f£o reveal that the peak splitting 
shows up at intermediate fields for &f£o ^ 5- Since &f£o 
is of the order of unity for the high-T c materials, the 
above mechanism may well explain the tunneling spec- 
trum observed at H = 6T . A detailed experiment on 
the field dependence may be worth carrying out. 

We finally consider the mixing of a d xy -w&ve or an s- 
wave component into the dominant d x 2_ y 2 wave to study 
how it affects physical properties in finite magnetic fields. 
Such a state has been suggested to be realized in finite 
fields at low temperatures accompanying a phase transi- 
tion [ p~9|j20| . Franz and Tesanovic treated this problem in 
the zero-field limit of an isolated vortex, where parame- 
ters were chosen in such a way that d x 2_ y 2+id xy state is 
already formed in zero magnetic field. Clearly this is not 
consistent with various experiments J^-|J] which indicate 
a d x 2_ y 2 pairing in zero field. We hence choose the pa- 
rameter g xy or g s in Eq. (Q) in such a way that a pure 
d x 2_ y 2 state is stabilized in zero field. More specifically, 

T cO V) / T c (T$/T c ) is put 0.31 (0.36) in the following d xy - 
wave (s-wave) study, whereas the critical value is 0.33 
(0.38) above which d x 2_ y 2+id xy (d x 2_ y 2 + is) state is 
formed in zero field below a certain critical temperature 
T (xy) ^ T (s)^ lower than T (xy) ( T (s)y Notice incidentally 



that T r = T r 



(x 2 -y 2 ) 



Figure 3 shows the field dependence 
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FIG. 3. Field dependence of the ratio | A xy \/\ A x 2_ y 2 \ 
(|A s |/|A :c 2_ !/ 2 1) at the center of the vortex-core square with 
T/T c = 0.1, fc F £ = 5, and T% y) /T c = 0.31 (T c ( s) /T c = 0.36). 
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FIG. 4. The tunneling conductance at the core center for 
the d x 2_ y 2+ id xy state (solid line) compared with that of the 
pure d x 2_ y 2 state (dashed line). Parameters are B / H C 2 = 0.3, 
T/T c = 0.1, fc P £ = 5, and T% v) /T c = 0.31. 



of the maximum-gap ratios |A a ; y |/|A :I .2_ y 2| and |A S | 
/| A. a .2_ y 2 1 at the midpoint of the next-nearest-neighbor 
vortices of the square lattice. One can see clearly that the 
component IA^jJ develops gradually without any transi- 
tions with the field strength. This fact can also be re- 
alized on a purely group-theoretical ground as follows. 
There are four mirror reflections in zero field which dis- 
tinguish the d xy state from the d x 2_ y 2 state. Those 
symmetry operations, however, no longer exist at finite 
fields due to the supercurrent flowing in the system. It 
hence follows that the mixing starts immediately above 
H c i. The same is true for the s-wave mixing as already 
pointed out p5| , p6| . The result implies that the reason 
for the observed plateau in thermal conductivity Jl9],|lJ , 
if it does exist conclusively, should be sought elsewhere. 
Finally, Fig. 4 compares the tunneling conductance of 
the d x 2_ y 2 + id xy state (solid line) with that of the pure 
d x 2_ y 2 state (dashed line) at T/T c = 0.1. We see lit- 
tle difference between the two even for B / iJ C 2 = 0.3 and 
^fCo = 5. A similar result has been found for the case 
of the s-wave component. These results suggest that the 
mixing without a transition brings only small changes for 
the physical properties in finite magnetic fields. 

In summary, we have carried out fully self-consistent 
calculations of the BdG equations for a basic d-wave 
model with continuously varying the field strength. We 
have thereby presented several new results connected 
closely with the detailed nature of the d-wave quasiparti- 
cles in the vortex lattice. The method used here will 
be helpful for the microscopic understanding of other 
anisotropic superconductors in finite magnetic fields as 
well as rotating Fermi superfluids. 

One of the author (T.K.) is indebted to Ch. Renner 
for a useful discussion. Numerical calculations were per- 
formed on an Origin 2000 in " Hierarchical matter analyz- 
ing system" at the Division of Physics, Graduate School 
of Science, Hokkaido University. 
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